clear matrix
clear mata
clear
set more off
set maxvar 10000

/*Data is sourced from Nigeria Demographic and Health Surveys https://www.dhsprogram.com/
Power Plants data is purchased https://www.africa-energy.com/live-data*/

global data "C:\Users\taiwo\Downloads\Nigeria Data\New_Analysis061424"
cd "$data"
use "$data\newdata_matching.dta", clear


global biophysical1 veg_index proximity_to_protected_areas proximity_to_water road slope


keep if treat_20 == 1
gen treatment = .
replace treatment=1 if distance <= 10000
replace treatment=0 if distance > 10000


gen treattime_10 = syear-pyear
replace treattime_10=. if treatment==0
replace treattime_10=. if treatment==.


*Used intervention group because these are the only househoulds who were treated

***Perform matching***
teffects nnmatch (wealthscore hhdedu hheadage hheadsex)  (intervention), ematch(dhsyear urban) osample(nomatch) ///
gen(stub2) atet biasadj(hhdedu hheadage) vce(robust)


***Generate Table A2***
tebalance sum

drop if stub21==.

gen id = _n

*Save only treat and matched control #2

save "$data\newmatchint.dta", replace



keep if intervention==1
keep stub21
rename stub21 id

sort id
quietly by id:  gen dup = cond(_N==1,0,_n)
keep if inlist(dup,0,1)
drop dup
save "$data\newctrlid.dta", replace

use "$data\newmatchint.dta",clear

keep if intervention ==0


merge 1:1 id using "$data\newctrlid.dta"
keep if _merge==3
drop _merge

save "$data\newctrlid.dta", replace
use "$data\newmatchint.dta",clear

keep if intervention==1 
append using "$data\newctrlid.dta"


save "$data\matched.dta", replace

use "$data\matched.dta", clear




teffects nnmatch (wealthscore hhdedu hheadage hheadsex)  (intervention), ematch(dhsyear urban) osample(nomatch2) ///
gen(stub5) atet biasadj(hhdedu hheadage) vce(robust)
***Rosenbaum Test
rbounds wealthscore, gamma(1.1(0.1)2)




balancetable treatment wealthscore hhdedu hhnumber hheadage hheadsex under5 urban $biophysical1  using "wealthm_102724.xls" if intervention == 0, vce(cluster projectid) replace

*eventdd
eststo binned: eventdd wealthscore hhdedu hhnumber hheadage hheadsex urban under5 $biophysical1  i.dhsyear i.state, timevar(treattime_10) cl(projectid) ///
accum leads(5) lags(4) method(ols) graph_op(ytitle("ATT") name(matched, replace) xtitle("periods away from power plants intervention") yscale(range(-0.5 1.1))) 
estat eventdd


***Generate Figure 2***
grc1leg matched, legendfrom(matched) 
addplot 1: ,  yscale(range(-0.6 1.0)) ylabel(-0.6(0.2)1.0) title("Wealth scores for different periods after matching", size(4) color(black)) legend(off) 
/*graph display, ysize(12) xsize(14)*/

***Generate Figure 3***

*csdid
gen year_csdid=0
replace year_csdid=pyear if treatment==1
csdid wealthscore hhdedu hhnumber hheadage hheadsex urban under5 $biophysical1 i.dhsyear i.state, time(syear) gvar(year_csdid) cluster(projectid) method(reg) notyet long2
estat all
estat event, window(-5 5) 
csdid_plot, name(Not_yet, replace) title("(a) Not yet and never treated")

csdid wealthscore hhdedu hhnumber hheadage hheadsex urban under5 $biophysical1 i.dhsyear i.state, time(syear) gvar(year_csdid) method(reg) cluster(projectid) long2
estat all
estat event, window(-5 5) 
csdid_plot, name(Never_treated, replace) title("(b) Never treated only")

*Combine plots
graph combine Not_yet Never_treated, name(combined_plot, replace) xsize(30) ysize(20)
graph display, ysize(14) xsize(20)


***Generate Figure 4***
csdid wealthscore hhdedu hhnumber hheadage hheadsex urban under5 $biophysical1 i.dhsyear i.state, time(syear) gvar(year_csdid) cluster(projectid) method(reg) estore(csdid3) notyet long2
estat all
estat event, window(-5 5) 
csdid_plot
csdid_estat event, window(-5 5) estore(csdid3)
estimates restore csdid3
local opts delta(rm) mvec(0(0.01)0.1) l_vec(l_vec)
local plot coefplot xtitle(M) ytitle(95% Robust CI)
matrix l_vec = 0\ 0 \ 0 \ 0 \ 1
honestdid, pre(3/7) post(8/12) `opts' `plot'
graph save honestdid_plot1, replace
*graph export coefplot.pdf, replace

csdid wealthscore hhdedu hhnumber hheadage hheadsex urban under5 $biophysical1 i.dhsyear i.state, time(syear) gvar(year_csdid) cluster(projectid) method(reg) estore(csdid6) long2
estat all
estat event, window(-5 5) 
csdid_plot
csdid_estat event, window(-5 5) estore(csdid6)
estimates restore csdid6
local opts delta(rm) mvec(0(0.01)0.1) l_vec(l_vec)
local plot coefplot xtitle(M) ytitle(95% Robust CI)
matrix l_vec = 0.5\ 0 \ 0 \ 0 \ 0.5
honestdid, pre(3/7) post(8/12) `opts' `plot'
graph save honestdid_plot2, replace

graph combine honestdid_plot1.gph honestdid_plot2.gph, title("Honest DID Analysis")
graph display, ysize(14) xsize(24)

sum installedcapacitymw, detail
g capacity_q4=(installedcapacitymw>r(p75))
g capacity_q3=(installedcapacitymw>r(p50))
g capacity_q2=(installedcapacitymw>r(p25))
g capacity_p90=(installedcapacitymw>r(p90))

_pctile installedcapacitymw, p(80)
gen capacity_p80 = (installedcapacitymw > r(r1))


eststo binned: eventdd wealthscore hhdedu hhnumber hheadage hheadsex urban under5 $biophysical1  i.dhsyear i.state if capacity_q3 == 0, timevar(treattime_10) cl(projectid) ///
accum leads(5) lags(4) method(ols) graph_op(ytitle("ATT") name(below_med, replace)  xtitle("periods away from power plants intervention") yscale(range(-0.5 1.1))) 
estat eventdd
eststo binned: eventdd wealthscore hhdedu hhnumber hheadage hheadsex urban under5 $biophysical1  i.dhsyear i.state if capacity_q3 == 1, timevar(treattime_10) cl(projectid) ///
accum leads(5) lags(4) method(ols) graph_op(ytitle("ATT") name(above_med, replace) xtitle("periods away from power plants intervention")  yscale(range(-0.5 1.1))) 
estat eventdd

***electricity

eststo electricity: eventdd wealthscore hv206 hhdedu hhnumber hheadage hheadsex urban under5 $biophysical1  i.dhsyear i.state, timevar(treattime_10) cl(projectid) ///
accum leads(6) lags(4) method(ols) graph_op(ytitle("wealth score") name(Binned_Years, replace) xtitle("periods away from power plants intervention") yscale(range(-0.5 1.1))) 
estat eventdd

eststo electricity: eventdd hv206 hhdedu hhnumber hheadage hheadsex urban under5 $biophysical1  i.dhsyear i.state, timevar(treattime_10) cl(projectid) ///
accum leads(6) lags(4) method(ols) graph_op(ytitle("Electricity Probability") name(electricity, replace) xtitle("periods away from power plants intervention") yscale(range(-0.5 1.1))) 
estat eventdd

***Generate Figure 5***
grc1leg below_med above_med electricity, legendfrom(above_med) 

addplot 1: , title("(a) Power plants below median capacity", size(4) color(black)) legend(off)
addplot 2: , title("(b) Power plants above median capacity", size(4) color(black)) legend(off) 
addplot 3: , yscale(range(-0.4 0.4)) ylabel(-0.4(0.2)0.4) title("(c) Electicity access likelihood", size(4) color(black)) legend(off)

graph display, ysize(14) xsize(20)

***Generate Figure A1***
***usual resident***

eststo binned: eventdd wealthscore hhdedu hhnumber hheadage hheadsex urban under5 $biophysical1  i.dhsyear i.state  if usual == 1, timevar(treattime_10) cl(projectid) ///
accum leads(6) lags(4) method(ols) graph_op(ytitle("ATT") name(usual_residents, replace) xtitle("periods away from power plants intervention") yscale(range(-0.5 1.1))) 
estat eventdd
grc1leg usual_residents, legendfrom(usual_residents) 
addplot 1: , yscale(range(-0.4 0.6)) ylabel(-0.4(0.2)0.6) title("Wealth scores for different periods using usual residents", size(4) color(black)) legend(off)


***Generate Figure A2***
***check electrical appliances ownership for households 
eststo radio: eventdd hv207 hheadedu hheadage hheadsex urban under5  $biophysical1 i.dhsyear i.state, timevar(treattime_10) ///
accum leads(6) lags(4) method(ols, cl(projectid)) graph_op(xtitle("periods before and after power plants")ytitle("probability of radio ownership") name(radio, replace) ytitle("radio") xlabel(-6(1)4))
estat eventdd
eststo television: eventdd hv208 hheadedu hheadage hheadsex urban under5 $biophysical1 i.dhsyear i.state, timevar(treattime_10) ///
accum leads(6) lags(4) method(ols, cl(projectid)) graph_op(xtitle("periods before and after power plants")ytitle("probability of television ownership") name(television, replace) ytitle("television") xlabel(-6(1)4))
estat eventdd
eststo refridgerator: eventdd hv209 hheadedu hheadage hheadsex urban under5 $biophysical1 i.dhsyear i.state, timevar(treattime_10) ///
accum leads(6) lags(4) method(ols, cl(projectid)) graph_op(xtitle("periods before and after power plants")ytitle("probability of refridgerator ownership") name(refridgerator, replace) ytitle("refridgerator") xlabel(-6(1)4))
estat eventdd

gen cooking_fuel = 1 if hv226 == 1
replace cooking_fuel = 0  if cooking_fuel == .
eststo cooking_fuel: eventdd cooking_fuel hheadedu hheadage hheadsex urban under5 $biophysical1 i.dhsyear i.state, timevar(treattime_10) ///
accum leads(6) lags(4) method(ols, cl(projectid)) graph_op(xtitle("periods before and after power plants")ytitle("probability of cooking fuel is electricity") name(cooking_fuel, replace) ytitle("electricity as cooking fuel") xlabel(-6(1)4))
estat eventdd

**Results no longer intuitive except for refridgerator
grc1leg radio refridgerator television cooking_fuel, legendfrom(radio) 

addplot 1: , yscale(range(-0.2 0.2)) ylabel(-0.2(0.1)0.2) title("(a) Radio", size(4) color(black)) legend(off)
addplot 2: , yscale(range(-0.1 0.2)) ylabel(-0.1(0.1)0.2) title("(c) Refridgerator", size(4) color(black)) legend(off)
addplot 3: , yscale(range(-0.2 0.3)) ylabel(-0.2(0.1)0.3) title("(b) Television", size(4) color(black)) legend(off) 
addplot 4: , yscale(range(-0.01 0.08)) ylabel(-0.01(0.05)0.08) title("(c) Electricity as Cooking Fuel", size(4) color(black)) legend(off) 

graph display, ysize(14) xsize(20)


use "$data\matched.dta", clear
sum installedcapacitymw, detail
keep  dhsclust dhsyear year  newyear pyear syear treatment treattime_10 state mean_temp veg_index travelminutes proximity_to_protected_areas proximity_to_water road slope projectid capacity_q3
duplicates drop dhsclust dhsyear, force
save "$data\matchedclust.dta", replace

use "$data\allyrwomen-var_int.dta", clear
append using "$data\allyrmen-var_int.dta"
merge m:1 dhsclust dhsyear using "$data\matchedclust.dta"
keep if _merge==3
drop _merge
global biophysical1 veg_index proximity_to_protected_areas proximity_to_water road slope
**encode sstate_1, generate(state)
recode v717 (6=1)(1/5=0)(7/11=0), gen(household_domestic)
gen some_education=primary
replace some_education=1 if secondary==1
gen manual=skilled_manual
replace manual=1 if unskilled_manual==1
gen sales_services=sales
replace sales_services=1 if services==1
replace sales_services=1 if household_domestic==1
replace sales_services=1 if clerical==1


***Generate Figure 6***

***manual
eventdd wealthscore hhnumber age education sex christian urban $biophysical1 i.dhsyear i.state if manual == 1, timevar(treattime_10) cl(projectid) ///
accum leads(6) lags(4) method(ols) graph_op(ytitle("ATT") name(manual, replace) xtitle("periods away from power plants intervention")  yscale(range(-0.5 1.1))) 
**sales_services
eventdd wealthscore hhnumber age education sex christian urban $biophysical1 i.dhsyear i.state if sales_services == 1, timevar(treattime_10) cl(projectid) ///
accum leads(6) lags(4) method(ols) graph_op(ytitle("ATT") name(sales, replace) xtitle("periods away from power plants intervention")  yscale(range(-0.5 1.1))) 
**agric
eventdd wealthscore hhnumber age education sex christian urban $biophysical1 i.dhsyear i.state if agric== 1, timevar(treattime_10) cl(projectid) ///
accum leads(6) lags(4) method(ols) graph_op(ytitle("ATT") name(agric, replace) xtitle("periods away from power plants intervention")  yscale(range(-0.5 1.1))) 

grc1leg manual sales agric, legendfrom(manual) 

addplot 1: , title("(a) Manual Labor", size(4) color(black)) legend(off)
addplot 2: , title("(b) Sales and Services", size(4) color(black)) legend(off) 
addplot 3: , title("(c) Agriculture", size(4) color(black)) legend(off) 
graph display, ysize(14) xsize(20)


